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Abstract 

Multiple measurement vector (MMV) problem addresses the recovery of a set of sparse signal vectors 
that share common non-zero support, and has emerged an important topics in compressed sensing. 
Even though the fundamental performance limit of recoverable sparsity level has been formally derived, 
conventional algorithms still exhibit significant performance gaps from the theoretical bound. The main 
contribution of this paper is a novel concatenate MMV and boost (CoMBo) algorithm that achieves 
the theoretical bound. More specifically, the algorithm concatenates MMV to a larger dimensional 
SMV problem and boosts it by multiplying random orthonormal matrices. Extensive simulation results 
demonstrate that CoMBo outperforms all existing methods and achieves the theoretical bound as the 
number of measurement vector increases. 

Index Terms 

Compressed sensing, multiple measurement vector, S-OMP, ReMBo, QR factorization, concatenate 
and boost 

I. Introduction 

Compressed sensing (CS) theory [1-3] addresses accurate recovery of unknown sparse signals from 
underdetermined linear measurements, and has become one of the main research topics in signal process- 
ing community. Most of the compressed sensing researches have been developed to address the following 
single measurement vector (SMV) problem [1-3]: 

(PO) : min |lx|[o, subject to y = Ax, (1) 

X 
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where y G M"*, A G R"*^'*, x G M", and ||x||o denotes the number of non-zero elements in the vector 
X. Since (PO) is a combinatorial optimization problem, alternative approaches such as greedy method [4, 
5] or basis pursuit [2, 6] has been widely used. Specifically, the basis pursuit addresses the following li 
minimization: 

(PI) : min ||x||i, subject to y = Ax, (2) 

where || • ||i denotes an li norm. The beauty of compressed sensing is that (PI) provides the exactly same 
solution as (PO) if the restricted isometry property (RIP) is satisfied [3]. More specifically, let k denote 
the sparsity of the solution vector x, i.e. ||x||o = k. Then, the mapping A satisfies the RIP if there exist 
a constant such that 

(l-efe)||v||2 < ||Av||2 < (l + efc)||v||2 , (3) 

for all nonzero vector v such that ||v||o < A;. It has been shown that for many class of random matrices, 
the RIP is satisfied with extremely high probabihty if the number of measurement is given by m > 
Ciklog{n/k) for some constant Ci [3]. 

Another important area of compressed sensing research is so called multiple measurement vector 
problem (MMV) [5,7-9]. The MMV problem addresses the recovery of a set of sparse signal vectors 
that share common non-zero support. Specifically, let 7^(X) denote the number of rows that have non-zero 
elements. Then, MMV problem addresses the following: 

min7^(X), subject to Y = AX, (4) 

X 

where Y G M"*^'', X G M"^'', and r denotes the number of measurement vectors. In MMV, thanks to 
the common sparse support, it is easily expected that the number of recoverable sparsity level 7?.(X) 
may increase with the increasing number of measurement vectors. This has been formally proven in the 
deterministic setup [5,7]. More specifically, Eq. (4) has the unique solution if 

7^(X) < ^{spark{A) + rank{Y) - 1) , (5) 

where spark{A) is the smallest possible number such that there exists a subgroup of spark{A) columns 
from A that are linearly dependent [10], and rank{Y) denotes the rank of Y that may increase with 
the number of measurement vectors. Note that Eq. (5) is a deterministic bound that guarantees 100% 

successful reconstruction. As a reconstruction algorithm for MMV problem, greedy algorithm such as 
S-OMP [5, 11], h convex relaxation method [12, 13], or M-FOCUSS [7] have been proposed. Recently, 
the novel algorithm called REduce MMV and BOost (ReMBo) has been proposed that outperforms the 
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conventional methods [8]. However, the performance improvement of the existing MMV algorithms are 
still limited and there still exist significant performance gaps from the deterministic bound in Eq. (5). 

The main contribution of this paper is to overcome the drawbacks of the existing method and to 
achieve the optimal performance. The new algorithm is based on the observation that QR factorization 
of the nonzero sparse matrix of X provides a sparse coefficient matrix. Then, by exploiting the sparsity 
of coefficient matrix, we can convert the MMV problem into a higher dimensional SMV problem, 
which can be solved using li minimization methods to estimate the common sparse support. Since a 
QR factorization of the solution matrix is unknown a priori, we boost the MMV problem by multiplying 
a set of random orthonormal matrices. As will be demonstrated later, the new algorithm can achieve the 
optimal performance as the number of measurement vectors increases. 

II. Conventional MMV Algorithms 

Simultaneous orthogonal matching pursuit (S-OMP) is a simple extension of OMP algorithm for MMV 
problem [5, 11]. More specifically, S-OMP chooses the maximum correlation between the residual and 
atoms: 

Algorithm S-OMP 

1) Set the residual matrix Rq = Y, the index set Aq = 0, and the iteration count t = 0. 

2) Find the index it that solves the optimization problem: 

max ||af^Rt||p (6) 

i 

where || ■ ||p denotes the Ip norm for 1 < p < 2. 

3) Set At = At_i U {it}. 

4) Determine the orthogonal projector Pt = Aa^ (A^ AaJ ^ A^^ onto the span of the atoms indexed 
in At. 

5) Calculate the new residual: 

Rt = (I-Pt)Y (7) 

6) Increase the step count t and go to step 2) if ||Rt||F > e- 

Even though S-OMP exploits the MMV structure, the performance dependency on the number of mea- 
surement vector has not been proven in the deterministic setup [5]. Recently, Gribonval et al [14], Eldar 
and Rauhut [15] prove the performance improvement using an average case analysis. 

In contrast, REduce MMV and Boost (ReMBo) addresses the MMV problem by reducing it to an 
SMV problem [8]. More specifically, MMV problem is first converted into an SMV problem by applying 
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a random vector w to the right hand side of Eq. (4) (reducing MMV). The resultant SMV problem is 
solved by high performance SMV algorithm such as basis pursuit. This procedure is repeated with other 
random vectors (boosting) until the termination condition is satisfied. 
Algorithm ReMBo 

1) Initialization: set iteration count t = and threshold value e. 

2) While (t < Maxlter) do 

a) Draw a unit norm random vector w G M''. 

b) Reduce MMV to SMV problem by multiplying w: 

y = Yw = AXw = Ax . (8) 

c) Solve Eq. (8) using SMV technique to get x. Let A denote the nonzero support of x. 

d) if (||x||o < k) and (|| y — Ax ||2< e) then stop, else t = t + 1. 

3) denotes the submatrix of X by collecting rows whose indices belong to A. Then, the final 
MMV solution is given by: 

(9) 

X* = 0, z ^ A. 

Using extensive numerical simulation, ReMBo has been shown to overperform other existing MMV 
algorithms [8]. The theory behind performance improvement of ReMBo has been recently suggested by 
van der Berg et al [9]. More specifically, based on the poly tope geometry, they showed that the recovery 
performance of SMV can be improved by limiting the sign of the unknown solution x to optimal patterns 
[9]. In ReMBo, the unknown solution x can be obtained from X by multiplying a random vector w, i.e. 
X = Xw. Hence, we have more chance to choose better sign patterns for x. The probabihty of choosing 
the optimal sign pattern should increase with the number of measurement vector. Hence, we can expect 
the performance improvement of ReMBo with increasing number of measurement vectors [9]. 

III. Concatenate MMV and Boost Algorithm 

In practise, the recovery performance of S-OMP and ReMbo are often saturated after a certain number 
of measurement vectors and there usually exists a significant performance gap from the deterministic 

bound, as illustrated in Fig. l(a)(b). Such performance saturation of ReMBo has been suggested by van 
der Berg et al [9], who observed the difficulty in visiting all desired sign patterns by drawing a finite 
number of random vectors w. In fact, some specific sign patterns are more dominantly observed [9]. 
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Therefore, rather than reducing the MMV into SMV, we are interested in directly exploiting the structure 
of matrix X for MMV recovery. 

Suppose the unknown matrix X G M"^'' with Tl(X.) = k assumes the following form: 

Q{n— fc)xr 

where K € W^'' denotes non-zero submatrix, and O € is (n — k)xr zero matrix, respectively. 

Now, consider a QR decomposition of matrix K. If r < k, the resultant QR decomposition is given by 



X 



(10) 



K = QR 



qi 



qr 



X 







X 



X X 



X 
X 
X 



(11) 



where Q 



qi 



denotes an r x r orthonormal matrix and R € M*^^ is a coefficients matrix 
from Gram-Schmidth orthonormalization, respectively. From Eq. (11), we can easily see that the number 
of non-zero elements in the coefficient matrix R is at most kr — r{r — l)/2. Second, if r > fc, the 
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X 





K = QR 



qi 



X 





(12) 











where the coefficient matrix R S W^^ has at most k{k + l)/2 non-zero elements in this case. 
In general, MMV problem can be represented as 

Y = AX = ASQ^ 

where Q is the r x r orthonormal matrix from QR decomposition, and sparse matrix S is given by 

Q(n— fc)xr 

Now, multiplying Q to the both sides of Eq. (13) provides us 



(13) 



(14) 



YQ = AS 



(15) 



Let Vec{Z} the vector of concatenated columns of a matrix Z. Using this, we can expand MMV problem 
to an SMV problem: 

YQ = AS ^ (r^'' A)Vec{S} = Vec{YQ} , (16) 

where F^'' denotes the r x r identity matrix. Furthermore, the structure of R matrix in Eq. (11) and 
Eq. (12) tell us that the total sparsity of Vec{S} is upper limited by: 

(fe-^), ifl<r<fe; 
iir>k. 

The same analysis can be applied even for more general situation where the location of non-zero rows 
are not clustered. Specifically, for general matrix X with 7^(X) = k, there always exist a permutation 
matrix P G M'^'''' such that 

Q(n— fc)xr 

Hence, the corresponding SMV problem can be easily obtained similar to Eq. (16): 



|VEC{S}||o < 



X = P 



PSQ^ . 



(18) 



(r^^ ® a)vec{ps} = vec{yq} 



(19) 
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Since the permutation does not change the sparsity level, we can easily see that the total sparsity of 
Vec{PS} is also upper limited by Eq. (17). Of course, the original MMV problem can be concatenated 
to a high dimension SMV without multiplying Q matrix as 

(F^'" ® A)Vec{X} = Vec{Y} . (20) 

However, in this case the total sparsity of Vec{X} is upper bounded by kr. Figs. 2(a)(b) illustrates 
the total and average sparsity versus the number of measurement vectors for the concatenated SMV 
problems Eq. (16) and Eq. (20). Average sparsity is defined as the total sparsity divided by the number 
of measurement vectors. We can easily see that the formulation Eq. (20) has no gain whereas the average 
sparsity keeps reduced with increasing r for the case of Eq. (16). 



(V^- ®A)VEC{X} = VEC{Y} 




(a) (b) 

Fig. 2. (a) Total sparsity and (b) average sparsity versus the number of measurement vectors when k = 10. 

However, the main problem of the proposed approach is that the orthonormal matrix Q is not known 
a priori since the computation of Q requires the exact knowledge of the unknown matrix X. Therefore, 
we randomly generate multiple candidates of orthonormal matrices and choose the one that yields 
best reconstruction performance. First, note that the sparsity level does not increase by multiplying Q. 
Furthermore, the resultant sparsity level should be in between that of Eq. (16) and Eq. (20) as illustrated 
in Fig. 2(a)(b), since random Q matrices approximates the orthonormal matrix in QR decomposition of 
unknown X. The pseudo-code implementation of the proposed algorithm is then as follows: 

Algorithm CoMBo 

1) Initialization: set iteration count t = and A = {1, 2, ... , n}. 

2) While (t < Maxlter) do 
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a) Draw a random orthonormal matrix Q G M''^^. 

b) Convert MMV to SMV problem by multiplying Q and concatenating the columns of YQ: 

Vec{YQ} = (F^'" A)Vec{S} (21) 

c) Solve Eq. (21) using SMV technique to get the estimate S. Let A denote the set of indices 
of the row vectors of S whose I2 norm is among the m largest. 

d) A ^ the set of non-zero row indices of (A^)^Y. 

e) If |A| < |A|, then A ^ A. Set i = i + 1. 
3) Let denotes the submatrix of X by collecting rows whose indices belong to A. Then, the final 

MMV solution is calculated as: 

X^ = At Y 

^ (22) 

Note the similarity of ReMBo and CoMBo pseudo-code implementation. In CoMBo, rather than reducing 
the MMV to SMV, we expand an MMV to a higher dimensional SMV problem by multiplying random 
orthonormal matrix Q and concatenates the resulting column vectors. Besides the slight differences in 
termination condition, the remaining steps of CoMBo are basically same as those of ReMBo. Hence, we 
conjecture that the role of random Q may have the same role to assign the best sign pattern. Specifically, 
CoMBo may exploit the possible combination of all sign patterns simultaneously by concatenating to a 
higher dimensional vector rather than reducing to the same dimensional vector. 

IV. Simulation Results 

In order to quantify the performance, we declare an algorithm is sucessful, if the following condition 

is satisfied: 

IIX-XII2 < 10-^||X||2, (23) 

where X denotes the ground-truth and X is its estimate, respectively. The simulation parameters are as 
following: m = 20, n = 30 and the number of measurement vectors are r = 1,2, 8,16, A; = 1 ~ 20, 
respectively. For ReMBo, we set e = 10~^. The Maxlter is set to 20 for ReMBo and 5 for CoMBo, 
respectively. The SMV algorithm used in both ReMBo and CoMBo is ^i-MAGIC [16]. All elements of 
sensing matrix A and the nonzero elements of X are drawn from the zero mean unit variance Gaussian 
distribution. Then, each column of A are normalized to have an unit norm, and the indices of nonzero 
values of X were chosen randomly. The success rate are averaged for 500 experiments. 
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Figs. l(a)(b) illustrate the recovery performance of S-OMP and ReMBo along with the deterministic 
bound in Eq. (5), respectively. Even though we can observe the performance improvement of S-OMP and 
ReMBo with increasing number of measurement vectors, there still exist significant performance gaps 
from the deterministic bound especially for large r. However, such performance gap can be significantly 
reduced using the CoMBo algorithm as shown in Fig. 3(a), which clearly show that the recovery 
performance of CoMBo algorithm approaches the bounds as r grows. Furthermore, with high probabilities, 
CoMBo even outperforms the deterministic bound of Eq. (5) for the case of r > 8. Similar average 
performance behavior has been observed in SMV problems [17], which shows that the sparsity more than 
spark{A) /2 can be still recovered with high probability. Fig. 3(b) illustrates that CoMBo outperforms 
other conventional MMV algorithms. Here, we also illustrate the performance of simple concatenation in 
Eq. (20). Thanks to the concatenation, the average performance is slightly better than ReMBo. However, 
the deterministic and average performance of CoMBo outperform these simple concatenation, which 
confirms the importance of boosting step. 



(FX'' ® A)VEC{X} = VBC{Y} 




(a) (b) 

Fig. 3. (a) Recovery rates of CoMBo, and (b) comparison of various MMV algorithms wtien r = 5. Ttie vertical lines of (a) 
denotes the deterministic upper limit of recoverable sparsity (Eq. (5)). 



V. Conclusion 

The conventional MMV algorithms such as S-OMP and ReMBo do not achieve the fundamental 
performance bound even with the increasing number of measurement vectors. In order to overcome 
the drawbacks, this paper proposed a novel CoMBo algorithm that concatenates MMV to a larger 



IEEE SIGNAL PROCESSING LETTERS, VOL. ?, NO. ?, JUNE 2009 



10 



dimensional SMV problem and boosts it by multiplying random orthonormal matrix. The CoMBo 
algorithm is a process of approximating the QR decomposition of the unknown solution matrix that 
allows reduction of average sparsity as the number of measurement vector increases. We compared the 
similarity between ReMBo and our CoMBo algorithm and showed that CoMBo is a generalization of 
ReMBo. Extensive simulation results demonstrate that CoMBo outperforms all existing methods and its 
performance approaches the deterministic limit as the measurement vector increases. 
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